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APPROXIMATE PREDICTIONS OF THE TRANSPORT OF THERMAL 

RADIATION THROUGH AN ABSORBING PLANE LAYER 

By Max. A. Heaslet and Frahklyn B. Fuller 

Ames Research Center 
Moffett Field, Calif. 

SUMMARY 

The transport of thermal radiation through an absorbing medium 
bounded by parallel walls is predicted . The walls are heated and are 
assumed to emit, reflect, and transmit the radiation isotropically. 

The medium is assumed to be in local thermodynamic equilibrium and to 
have radiative characteristics that can be averaged over the entire 
frequency range. First, relations are given whereby the solution 
becomes available once the problem associated with the case of opaque, 
black walls is solved. Second, different methods are employed to 
derive approximations of radiative flux and temperature distribution 
through the medium. Simple formulas, of interest in engineering design 
and analysis, appear as a by-product of the study. Comparisons are 
made with previously published numerical results. 

INTRODUCTION 

A number of references are now available in which the equations 
characterizing the transport of thermal radiation through an absorbing 
medium are formulated and mathematical techniques leading to actual pre- 
dictions are developed. Foremost among these are the treatises of 



Chandrasekhar [1], Kourganoff [2], and Sobolev [ 3 ]. The analysis 
retains considerable algebraic complexity, however, partially because 
of the multiplicity of physical parameters affecting the phenomena and 
the number of additional idealizations that are of practical interest. 
Thus, in spite of the theoretical understanding that has been achieved, 
and although numerical solutions by means of high-speed, computers offer 
no essential difficulties, the problem of efficient presentation of 
results remains a challenge. The present paper attempts to contribute 
some efficiency of presentation to predictions of radiative transport 
through an absorbing medium contained between two parallel walls with 
specified temperatures and known radiative characteristics. 

Two parallel objectives are to be kept in mind. First, a close 
connection is established between the general problem of radiative 
energy transport between walls that emit, transmit, and reflect iso- 
tropically and the problem of transport between opaque, black walls. 

The latter case has been treated by Usiskin and Sparrow [ 4 ] and numeri- 
cal solutions have been given for different values of the optical thick- 
ness of the plane layer between the walls. It will be shown here that 
the results of [4] may be applied with only minor modification to the 
more general case. Second, since the calculations of Usiskin and 
Sparrow are available for comparison, approximations to the solutions 
of the governing equations are studied. In this way, rather simple 
expressions are derived that provide a surprising accuracy and that 
certainly should be useful in preliminary design or engineering analysis. 


2 



In the next section the governing equations are derived. The more 
or less conventional approach leads to an expression for flux of radia- 
tion that was given explicitly "by R. and M. Goulard [5]. Differences 
in terminology were considered necessary and the derivation of wall con- 
ditions was cast in a different form, all of which prompted the inclu- 
sion of this basic material. The discussion of the general solution 
gives a proof that iterative methods can he used to solve the basic 
integral equation and then relates the problem to the case of black 
walls. The remaining portion of the paper deals with approximations 
and comparisons between solutions involving differing orders of accuracy. 

Three methods, yielding increasing degrees of accuracy, will evolve. 
The first, based on the use of an exponential influence function in 
place of the exact exponential -integral function, simplifies the analysis 
so much that all results can be expressed algebraically. The detailed 
distribution of the temperature (or the emission) within the medium suf- 
fers increasing inaccuracy near the walls but further integration to 
determine flux between the walls leads to predictions that are never in 
error more than a few percent. Perhaps more important is the fact that 
the formulas are easily manipulated to display major effects of varia- 
tions in the various physical parameters. The second method is similar 
in approach to the Milne -Eddington approximation used in the study of 
stellar atmospheres and improves the accuracy of the predictions without 
introducing “undue complexity. Finally, an iterative calculation of a 
particularly simple form provides estimates that agree well with the 
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more accurate numerical solutions that are available. These iterative 


calculations are then used to establish the degree of confidence one 
can have in the other approximations. 


TABLE OF SYMBOLS 


A 

En(x) 

H(S,| l ) 

l(x,d) 

<3 

J 

Ko, K x 
L 

m, n 


4 

r 

t 

T 

x 

Z 

a 


(1/2) [(1/2) - Es(S l )] 


Planck's function (see eq. (5)) 


integroexponential function of order n; En(x) 



-xt t -n 


dt 


l 

integral used in iterative solution (see eq. ( J+8b ) ) 
specific intensity, energy per unit area, time, solid angle 
local emission coefficient, per unit mass 
source function (see eq. (5)) 

weighted integrals of emission function (see eqs. ( 27 ) and ( 32 b)) 
geometric thickness of plane layer (fig. 1) 
coefficients used in exponential approximation of Es(x) 

(see eq. (35)) 

rate of energy transport per unit area 

reflectivity 

tr ansmi s s ivity 

temperature, absolute 

geometric depth in absorbing layer 

see equations ( 27 ) 

absorptivity 
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P(|) emission function (see eq. (ll)) 

7 coefficient of cubic term in approximation of emission 

function (see eq. ( 46 )) 

A see equation ( 29 ) 

e emissivity 

0 angle between ray and element normal to surface (see fig. 1) 

k local absorption coeff icient, per unit mass 

A slope of^linear approximation to emission function (see 

eq.. (45a)) 

H cos 0 

| optical depth in absorbing layer (d£ = pK dx) 

p local density of absorbing medium 

a Stefan 1 s constant 

ep(|) dimensionless form of emission function (see eq. (25)) 

Subscripts 

o evaluated at f* = 0 

r L 

L evaluated at I = = / P K 3 x 

Jo 

Superscripts 

+ right -going (| increasing) quantity 

left-going (£ decreasing) quantity 
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GOVERNING EQUATIONS 


Figure 1 shows the orientation of the coordinate system to he 
used. The two walls are situated at x = 0 and x = L. Let conditions 
at either wall he indicated hy the subscript i where i = 0 or L. 
Then, we have temperature, Tj_, emissivity, e±, absorptivity, a,j_, 
reflectivity, r-j_, and transmissivity, t-j_, prescribed as known boundary 
conditions where 


a i = € i 


(la) 

+ ti + r i = 1 


(lb) 


It is assumed that the gas and walls have characteristics that may be 
averaged over the entire frequency range of the radiation, so that a 
gray analysis applies, and the walls emit and reflect diffusely. 

The basic equation of radiative transfer is (see, e.g., [2]) 

P = -pKl(x,p) + pj (2) 

dx 

Equation (2) is preferably expressed in terms of the independent vari- 
able | (optical depth) where 

d| = P K dx (3) 

and thus becomes 

P — + I = J (4) 

at 
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vlier e J (=j/p^) is called the source function. Provided the refrac- 
tive index of the medium is 1 and local thermodynamic equilibrium 
is achieved, the source function is given by 


where 


r 00 

J = i *> 4v -i T ‘ 


By 


2h v 3 exp( -hv/kT) N 
c 2 1 - exp( -hv/kT)/ 


Planck’ s function 


(5) 


v frequency 

h,k, cr Planck, Boltzmann, and Stefan constants 

c velocity of light 

The notation 


i = i + (i,n) o < n < l 


(6) 


serves to distinguish the rate of flow of energy per unit area and 
solid angle for positive and negative values of p(=cos 0). With this 
notation, the solution of equation (4) in the two regimes becomes 


i + (g,n) 


l + (0)exp( -g/n) + 



J(l 1 )exp[-(| - g 1 )/ [ a] (dg^/p.) 


(7a) 


I (I,n) = l"(lL) ex P(-(^ - I L )/nl 



j(l 1 )exp[-(^ - | 1 )/n](d| 1 /-n) 

(7b) 


where I + (0) and I~(£l) are the isotropic specific intensities induced 

r L 

by the walls at g = 0 and g^ = pK dx, respectively. 
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The total rates of energy transport per unit area in the positive 
and negative x directions are 


q + (|) = 2rt f I + (|,p)p dp 

(8a) 

'-'o 




<T(i) = 2jt / I~(t,p)l-i dp 

(8b) 


The values of q. + (|) and q~(i) are sometimes referred to as the right 
and left fluxes of radiation. Equations (j) and (8) yield the two 
relations 


q + (l) = 2qJ EfeU) + 2rt T 1 J(l 1 )E 2 (| - ijd^ 

r fi L 

q-(l) = 2q£ E 3 (Z l - I) + 2rt / J(l 1 )E 2 (| 1 - gjd^ 


(9a) 

(9b) 


where q^ = jtI + (0) q£ = jrl”(^ L ) and E n (|) is the nth exponential 
integral defined as 

Ms) 


If local unidimensional energy flux (net energy transport per unit 
time and area) is denoted q(|), one has 

cl(S) = q + (l) - q-(i) (10) 

and q(|) >0 obviously corresponds to a net energy flow in the positive 
x direction. From equations (9) and (10) follows the fundamental rela- 
tion. After setting 
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itJ(i) = P(6) = aT 4 U) (11) 

one gets 

<l(l) = 2qJ Ea(l) - 2q.£ E 3 (l L - |) 

r lL 

+ 2 / p(g )sgn(g - Ii)E 2 (U - l-iDaij. (12) 

It will be noted that the expression sgn(g - | 1 )E 2 (I| - g-J)* in the 
integrand of (12), is discontinuous when g = g. This requires care 
in differentiation. 

Equation (12) is the flux equation for our unidimensional problem. 
The notation does not agree precisely with conventional astrophysical 
terminology but the choice of variables should preclude any ambiguities. 
The value of q(g) differs by a factor it from the astrophysicists 1 
flux, usually denoted E. Sign conventions are also changed. The inde- 
pendent variable g is, of course, dimensionless since the volumetric 
absorption coefficient pK is measured in terms of the reciprocal of 
radiation mean free path length. 

The dimensionless parameter g^ (optical thickness of layer) will 
be seen later to play a unique role in fixing the nature of the varia- 
tion of the emission function P(g). It should be remarked that in some 
references the volumetric absorption coefficient is assumed a constant 
so that g-Q becomes pftL. The transformation to optical path length 
in equation (3) shows that this restriction is not necessary. In a 
recent paper, Probstein [6] has developed an analogy between pkL and 
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the inverse Knudsen number of low-density fluid mechanics. The Knudsen 
number appears in molecular transport phenomena and is the ratio of 
molecular collision mean free path to the characteristic flow length. 
The term pKL can similarly be expressed as the ratio between charac- 
teristic length and the mean free path of a photon. This concept is of 
considerable heuristic value and enables Probstein to present, for the 
case of black walls, an approximate expression for flux by means of the 
analogy. Similar approximate results that apply to more arbitrary wall 
conditions will be considered in the final section of this paper. 

Calculation of Input Conditions at Walls 

The magnitudes of qj and q£ are related to the wall temperatures 
(or direct wall emission) as well as the wall parameters introduced in 
equations (l). Two equalities exist which serve to fix these values. 
Thus, one can write the two balances 

qj = € 0 oT 0 4 + r Q q£ (13a) 

= € L aT L 4 + r L^L (l3b) 

These relations simply state that the energy per unit time and area 
coming from each of the walls is equal to the direct emission plus the 
reflected portion of the incoming radiation from external sources, that 
is, from the other wall and the medium between the walls. The quantita- 
tive expressions for q£ and q£ are, from equations (9), 


I 


10 



I 




q£ = 2q£ E 3 (| l ) + 2 f PdjEaCl^d^ 

J o 


r lL 

c£ = 2qJ E 3 (l L ) + 2 / 3(i 1 )E 2 (| L - 

From equations (13) the desired conditions at the walls now follow; 
thus 

q+ = [1 . ltf 0 r L E 3 2 (£ L )] 1 € 0 aT 0 4 + 2r 0 e L E 3 (| L )oT L 4 

p^L p^L 

+ 2r 0 / PUjEsUJdi! + 4r 0 r L E 3 (| L ) / PUJE^Il - ijd^ 
J Q J 

(14a) 

l£ = [1 - 4r 0 r-[-E 3 2 ( 1^) ] 2r^€ Q E 3 ( |^)0T O 4 + e^cfT^ 4 

p 

P(i 1 )E 2 (i L - i 1 )d| 1 + 4r 0 r L E 3 (| L ) / PCljEsU^d^ 

(14b) 


+ 2r T 


II 


Basic Integral Equation for Constant Flux 

For constant wall temperatures and local thermodynamic equilibrium, 
q( | ) , as given in equation ( 12) , must he a constant. The derivative of 
q(|), from equation (12), is 


2^ = _ 2 qJ E 2 (|) - 2q£ E 2 (l L - |) + 4P(|) - 2 J Pd^E^ll - | 1 l)dg 1 

(15) 
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and for constant flux one gets the basic integral equation 


P(S) = \ <a+ e 2 (i) + | <a£ e £ (s l - 6) + | J P(s x )Ei(l ^ 

(16) 

Two special cases of equation (l6) are of particular interest. 
Firsts for opaque, Hack walls, tq =0, rq =0, and the wall emissions, 
as given in equations ( 14) , are 

q£ = 0T o 4 , q£ = oT L 4 

Usiskin and Sparrow [4] have given numerical results corresponding to 
these conditions. The importance of this case is enhanced once one 
observes in equations (l4) that q^ and q£ have no explicit dependence 
on | ■ and that the analysis can proceed formally without a precise 
stipulation concerning physical conditions at the walls. Equation (l6) 
is thus used to determine solutions in terms of the parameters q^ 
and q£ and the connection between these parameters and the actual 
boundary conditions is established later in an auxiliary calculation. 

A second case of interest arises in relating equation (l 6) to the 
problem of a plane, parallel, semi -infinite stellar atmosphere. This 
idealization occurs when the left wall at x = 0 is transparent (t 0 = 1, 
r 0 = 0) and the right wall is allowed to recede to an infinite distance. 
Here, the driving terms in the right member of equation (l 6) vanish, 
since qj = 0 and {-p, = oo. The integral equation becomes homogeneous 
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( 17 ) 


p(l) = \ J PdjEidl - 

and the physical constants that should determine the magnitude of the 
emission function 0(|) seemingly disappear from the problem. Actually, 
the one constant that fixes the radiation level in this case is the con- 
stant value of flux. The function 0(g) must he related to q which 
means that equation (12) cannot he disregarded. Equation (12) is, in 
fact, the fundamental governing equation; the derivation leading to 
equation (l 6), or ( 17 ), involves the loss of a physical constant, and 
it might he expected that under the most general conditions a restoration 
of the constant would he required. 

In astrophysics the study of a semi-infinite atmosphere under con- 
ditions involving thermodynamic equilibrium is referred to as Milne* s 
restricted problem, and equation (17) is called the Milne (or Milne - 
Schwarzs child) integral equation. A more than substantial literature 
exists on the mathematical analysis of this particular idealization. 
Milne* s equation is obviously a limiting case (g-^ -> 00 , q^ 0) of equa- 
tion (l6) and in this sense the highly accurate solutions in existence 
for the Milne problem can be used as a check on extreme conditions for 
the present study of radiation transport between parallel walls. 

Passage to the limit is not always direct, however. 
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GENERAL SOLUTION 


Equation (l6) is a Fredholm integral equation of the second kind. 
It takes the classical form 

r lL 

3(1) = f(s) + \ J 3(& 1 )E 1 (|g - IJJdlj, ( 18 ) 

vhere 

2f(l) = qj E £ (g) + q£ E 2 (| L - g) 

Mach of the mathematical literature devoted to this type of problem is 
built around the use of the Schwarz inequality (see, e.g., Tricomi [ 7 ], 
p. 50 et seq. ) which in turn leads to integration of the square of the 
kernel function. Rather than relate the analysis to general concepts, 
it is more expedient to derive what results are needed. The knowledge 
that the kernel Ei( | g - £ x | ) is everywhere positive will be used. 

Let g-^ be a finite constant. We prove first that the solution 
can be calculated by repeated iteration. The first (n + 1) iterations 
in the calculation of 3(g) lead to the sequence 3n+i(i) where 

3i(|) = f(S) 

Pa(6) =f(i) +| J f(l 1 )E 1 (lg - | 1 l)dg 1 (19) 


1^ 



p Il 

Pn+i(6) = f U) + | J fdjEiCll - SjHlj. + . . . 

1 r^ L p 

+ ^ / ai n Ei(l6 - i n l)/ ain.^idin 

2 'Jq ^0 

n^L 

/ d| 1 f(6 1 )B x (l6 2 - sj) 


fin-il) 


The relations 


f(6) < f. 


max 


(20a) 




Bids - liDagi 


-max 


1 

2 


2 - Efed) - E 2 U l - |) 


< M < 1 


-max 


(20b) 


apply for finite Then since the kernel is positive, repeated use 

of the mean value theorem for integrals yields 


P n+ i(l) < wd + M + M 2 + . . . + rf 1 ) < 

The sequence for positive terms therefore converges for all g. 

For a nonunique solution to exist, a finite, nonvanishing solution 
of 

P(s) = \ f L Pd^Eidi - liDdli (21) 

would be required. Meghreblian [8] has shown that this is impossible 
as follows. From equations (20b) and (21), the inequality 

1 P ^ 

lP(l)l <2 J I P(l x ) lEi(l 5 - I x I )d| x < M| Pmax I < I Pmax I ( 22 ) 

holds. But this cannot apply for all g in the range 0 < g < g^. 
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The finiteness of is an important condition underlying the 

above results. Retaining this condition and returning to equation (l6) 
one can show that 


p (v) ■ I K + 


(23) 


and that the function is antisymmetric about the point | = ljj/2. In 
fact, straightforward manipulation of equation (l6) reveals that the 
function P(|) + " 0 " + q£) satisfies the homogeneous form 

of the integral equation (eq. (21)). But for finite we have noted 

that the solution must vanish throughout the range of integration and 


this establishes the antisymmetry property 

qj + qf qj + qf 

P(6) l-S) 


(24) 


When | equation (23) is seen to hold. 

In actual calculation and in the presentation of results it is 
convenient to introduce the function cp(|), where 


<p(S) = 


ED 

% ~ 4 


+ *£ 
2 ^£ - O 


(25) 


From equations (23) and ( 2k ) one concludes that cp(|-^/2) = 0 and that 
cp(|) is antisymmetric about the point £ = 1^/2. The integral equa- 
tion for cp(|) is, from equations (l6) and ( 25 ), 

i i r ih 

cp(i) = i [E 2 (i L - |) - Efe(|)] + | J cp(6 1 )E 1 (| £ - 6 1 l)d| 1 ( 26 ) 

The function cp(|) has no explicit dependence on q^ and q£. 
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Calculation of Solution 


The main objective of the analysis was to calculate the emission 
function P(£), knowing the wall temperatures and radiative character- 
istics. Equation ( 25 ) shows that when P(|) is made dimensionless by 
dividing by q£ - q£, the problem separates into two calculations: 

(a) Determination of the universal functions 9(1) which depend 
on the single parameter 

(b) Determination of (q£ + q£) and (q£ - qj) which are constant 
for given conditions but which depend on all of the given 
parameters of the problem as well as the solution 9(1). 

Before the numerical solutions are presented, additional relations 
will be given to establish the interdependence between the various 
parameters. We assume here that has been fixed and that the 

function 9(1) has been calculated. Substitution from equation (25) 
into equations (l4) then yields, after setting 


Z - 1 - 4r G rj j E3 2 ( 1^) 




Ko(l L ) ^ 



cp(l 1 )E 2 (| 1 )d| 1 


( 27 ) 


the two simultaneous equations for 


<£ q£ 
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(28a) 


q£ [Z + 2r 0 (Ko - A) - 4r 0 r L E 3 ( g L ) (Kq + A) ] + q£ [-2r 0 (Ko + A) 

+ l4-r 0 r L E 3 (i L )(K 0 - A)] = e 0 oT 0 4 + 2r Q E 3 ( | L ) s L crT L 4 

qj [-2r L (Ko + A) + 4r 0 r L E 3 (| L ) (Kq -A)] + q£ [Z + 2r L (Ko - A) 

- ltf 0 r L E 3 (| L )(Ko + A)] = 2r L E 3 (^ L )e 0 aT 0 4 + 6 L oT L 4 (28b) 

If the additional notation 

A = Z + 2[(r 0 + rj j )(K 0 - A) - 4r 0 r^E 3 ( 1^) (Kq + A) - 8r 0 r£KoA] 

= 1 - r 0 r L + (r 0 + r L - 2r 0 r L ) 2Ko(g L ) + E 3 (| L ) - |J (29) 

is introduced, the determinant of the coefficients of the simultaneous 
equations is Z-A and the solution of equations (28) yields 

qj = ||e 0 tJT 0 4 [l + 2r L (Ko - A) ] + 2r 0 e L aT L 4 [E 3 ( | L ) + Kq + A]| (30a) 

q£ = i|2r Leo aTo 4 tE3(i L ) + K 0 + A] + € L oT L 4 [l + 2r 0 (Ko - A)]} (30b) 

Tims 

<l£ “ ^ = i [( ' 1 ~ r o) € L aT L 4 " ^ ' r L) € oC«o 4 J (31a) 

q£ + <3+ = ^{s L 0T l 4 [(1 + r 0 )'+ 4r 0 (K 0 - A)] 

+ € 0 ctT 0 4 [(1 + r L ) + 4r L (Ko - A)]j- (31b) 

From equations (12) and (25), at £ = 0 

1 = (q£ - q£) - Es(g L ) - 2 Ko(| l ) (32a) 
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and at | = 1^/2 

<1 = (q£ - q+) - 2 E 3 (~) + ^KxU l ) (32b) 

where 

n ^l/ 2 /|t \ 

Ki(6 L ) = J o (j- - i ly ) ai x 

Flux, q, and the emission function, p(|), are finally expressed 
in terms of the universal function cp(| ) , the boundary conditions involv- 
ing crT 0 4 , crT^ 4 and the wall parameters as follows: 

P(£) = ^{[(l “ r o)e L aT L 4 " ( x " r L )e 0 aT 0 4 ]cp(i) + | e L aT L 4 [(l + r Q ) 

+ ^r 0 (Ko - A)] +| e 0 tJT 0 4 [ (l + r L ) + 4r L (K 0 - A)]j- 

= ^ C E (1 - r o) s L (jT L 4 - (l-r L )e o 0T 0 4 M!) + | 6 L aT L 4 {l + 2r 0 [2K 0 (| L ) 

+ E 3 (i L )]} + | eo^o 4 ! 1 + 2r L [2K 0 (i L ) + E 3 (l L )]}) (33a) 

QL = ^ [1* + 2E3(|^) + ^K 0 (^l)][( 1 - r o) € L cr ^L 4 " ^ ~ r L^ e o°^o 4 ] 

= - | E S 00 - 2K ± [(1 - r 0 )e L aT L 4 - (l - r L )e 0 oT 0 4 ] (33b) 

where A, K^, Ki, and A are given in equations (27), (29), and (32). 

The integral equation (26) for cp(|) can be rewritten in an 
alternative form. From equations (12) and (32) this is 
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( 34 ) 


~ 2 ~ Es ^l) + ^ 3 ( 1 ) + E 3 (| L - l) 


p^L 

= 2 / cpCl^CsgnCl - SjEaCll - | x l) + E 2 (| 1 )]d^ 1 


' J o 

Equation (2 6 ) can "be rederived from equation ( 34 ) "by taking the deriva- 
tive of the latter with respect to | • Either equation is adaptable to 
the calculation of cp(g) by numerical means. The first is a Fredholm 
equation of the second kind and, as shown previously, the proof of the 
convergence of iterative methods offers no analytic difficulty. Numeri- 
cally, however, the singularity in the kernel at | is an incon- 

venience requiring special attention. Equation (34) is a Fredholm equa- 
tion of the first kind, is less well adapted to general analysis, but 
is numerically more tractable since the kernel is finite everywhere 
although it does possess a step discontinuity at g = | 1 . 

Equation (33a) shows that the variation of p(g) depends essentially 
on the variation in cp(g). For theoretical predictions, therefore, the 
basic calculations can be carried out for the case of opaque, black walls 
and for a given optical thickness g-^ in order to fix cp(g). Thus, the 
graphical results given by Usiskin and Sparrow [4] or by Meghreblian [8] 
can be used directly. The terminology used here differs from these ref- 
erences; it suffices merely to note, however, that cp(g) is the devia- 
tion of the dimensionless emission function from its value optically mid- 
way between the walls. Following this determination, the exact level 
of the function 0(g) can be found after an additional integration to 
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It -will be shown in the next section 


determine K 0 or K 0 and K x . 
that good approximations of the solutions can actually be carried out 
analytically. Graphical results will be given then. 

APmoxmtmi SOLUTIONS 

In this section, approximate solutions to the transfer equation 
are given in analytic form. First, the nature of cp(|) will be deter- 
mined when the kernel is arbitrarily replaced by an exponential function. 
This type of approximation is often used in more complex problems involv- 
ing the coupling of radiative and other modes of energy transport (see, 
e.g., [9])* Situations like the present one, in which different orders 
of accuracy can be achieved, are of value in assessing the degree of 
accuracy that can be expected. Second, the exact kernel is retained 
but the analysis is restricted to the determination of a linearized form 
of cp( g ) together with its improvement by means of a single iteration. 

In general, this solution represents an improvement in accuracy over the 
first case but at some cost in ease of calculation. Third, a more 
refined process, representing a modified iteration of the exact equation, 
will be carried out. These latter results compare well with the pre- 
viously published numerical solutions for 0 < ^ < 10. 

A commonly used approximation for E^Ce) involves its replacement 
by an exponential function. More specifically, the relation 

E 2 (i) » rne" n ^ (35) 
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is introduced where m and n are chosen in some plausible manner. A 
plot of the two functions shows that small deviations can he maintained 
throughout the full range of | , even though local discrepancies in the 
derivatives of the functions become excessive. This approximation was 
implicit in the early work of Eddington [10] on stellar atmospheres, 
the rationalization stemming from an averaging over the angular variable 
Q (or \i ) of the integrals appearing in the derivation of the transport 
equations. Different values of m and n are found in the literature 
since compromises must be made when mean values of specific intensity, 
flux, and radiation pressure are required. A discussion from the point 
of view of the astrophysicist is given by Ambartsumian ([11], p. 17 et 
seq.. ) . More recently, Vincent! and Baldwin [12] have proposed m = n 2 /3,> 
n =3 I.562 which follow after the requirement is imposed that flux be 
given properly in the Rosseland limit of strong absorption and that equa- 
tion (35) correspond to a least squares fit between the two functions. 

The first of these conditions is well founded physically and is common 
to all approximations. The second condition is obviously more arbitrary. 
Proposed values of n in the literature fall usually in the range 
3/2 < n < ; the value 3/2, preferred by most authors, is the one 

associated with Eddington 1 s approximation. 

If the flux equation is written as 

z ^ x = -CE 3 (0 + E 3 (£ l - I)] + 2 / qpClJsgnU - l 1 )E 2 (U - 
% " J ° 

( 36 ) 
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the exponential approximation can he introduced in the integrand from 
equation (35)- The resulting integral equation can he solved explic- 
itly hut the degree of approximation is affected only slightly and the 
analysis is made much simpler if the further approximation 

n °o 

*3(1) -f fcdjat, i = - | •"* (37) 

is introduced. The integral equation then becomes 


m 




e -n ^ + e - n(l L^) 


>^L 


2m f cp(| 1 )sgn(| - | 1 )e _n * d£ x 
d 0 

(38) 


For constant flux, the solution is 


<p(6) = - n ■(& - ^ 

2 + n| L 


(39) 


and consistent with the approximation used, one has 



- ,-sl) , 

,2 + n| L J 



2 

2 + n£ L 




For all values of the parameters the emission function P( | ) is therefore 
a linear function of | . From equations (31a) and (32b) 


-q. 4-n 

% ~ ' 3(2 + n| L ) 


(40a) 



1 + (r Q + r L ) 


(l - r 0 )e L aT L 4 - (l - r L )e 0 aT 0 4 

2[n - (3/2)] - n£ L [n + (3/2)31 2r r n 2 “ n ^L 
3(2 + n| L ) J r ° rL 32+ n| L 


(40h) 
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At extreme values of 


1-^, flux is given by 


-q. 




§L « 1 

i L >:> 1 


(41) 


When n = 3/2 and when the walls are opaque and black, equation ( 40 a) 
is the interpolation formula for flux derived hy Prohstein hy completely 
different methods. Prohstein* s graphical comparison of his formula and 
more accurate results indicated good agreement. Further comparisons, 
including different emissivities and reflectivities, will be given in 
figure 2. 

Equations (41) give intuitively obvious results: for optically 

thin media « l) flux is merely the difference between the total 
emissions from the two walls; for optically thick media, where the mean 
free path of the radiation is small, the flux is in agreement with the 
predictions of Rosseland* s approximation. The Rosseland theory reduces 
the radiation transport equation to the heat conduction equation with 
an effective, nonlinear conductivity. For a gray material with unit 
index of refraction, flux is given by 


q * 


4 

3 dt, 


(42) 


Equation (42) holds only for |t » 1 and boundary conditions are strictly 
applicable only to the case of opaque, black walls. Integration of equa- 
tion (42) for constant q then leads directly to the second of equa- 
tions (41) . 
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An improved approximation can be achieved through a study of equa- 
tion (26) directly. The function cp(|) is known to he asymmetric about 
| = |p/2. To a first order , therefore, it is proper to set 
cp(|) =7\[| - (|p/2)] and seek ?\. Direct substitution may then be made 
into equation ( 2 6 ), or into the derivative of this relation, which has 
the form 


q>'(S) = ^ [1 + 2cp(0)]Ei(£) + i [1 - ap(i L )]Bi(6 L - |) 


+ 




cp’UjEidi - i 1 l)ae 1 


(43) 


The calculation gives, at g = 1-^/2, 


/SlN , \ fih\ El ^l/ 2 ) 

[—■) - cp(0)E 1 


' (44) 


The linear approximation for cp yields cp(0) = This latter 

expression, together with equation (MO, gives the following predictions 


A = 


e SL/2 E 1 (i L /2) 


(45a) 


Cp(l) = 


i lL ^ 2 Ei(l L /2)[g- (| l /2)] 


Kq = A^[E 4 (0) - E 4 (g L )] - ^ [Ba(0) + E 3 (g L )] 


Ki = A ^ E 3 ) - 


E 4 (0) - E 4 ( ^ 


(45b) 


(45c) 
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_ a k r E 2 U l /2)E 3 (! l /2) - E 1 (| l /2)E 4 (| l /2) 

% ~ % 3 L EiUj/2) 

q£ - qj = [(1 - r 0 )€ L aT L 4 - (l - r L )e 0 aT 0 4 ] -jl 

+ (r Q + r L ) - \ - * (t " |) + C 1 - " 2 ^(l L ) 

-2r Q r L -7\ (^ ' f) + (l “ 7 s| Ii )E 3 ( i L ) - 2AE 4 (i L ) } (^5e) 

Equation (45a) is in agreement with the determination of slope given by 
Meghreblian [8] for black walls but equation (45d) does not conform to 
his result. The linear approximation of cp(|) does not yield constant 
flux and a decision must be made concerning its evaluation. The flux 
at as given in equation (45d), gives much better accuracy than 

flux evaluated at the wall. 

Equations (45b) and (45d) represent improvements over the simplifi- 
cations inherent in equations (39) and (40). The degree of accuracy of 
equation (45b) is, of course, limited since the predicted emission is 
expressed as a linear function of | throughout its range. By virtue 
of the additional integration to get flux, however, equation (45&) should 
provide good estimates. One notes that for extremes in |^, equa- 
tion (45d) reduces to the predictions of equations (4l), as indeed it 
should from physical considerations. 

An improved expression for cp(£) now follows by iteration when 
equation (45b) is used to evaluate the integral in equation (26). This 
process is much to be preferred, analytically, to the more general 
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iterative scheme employed in equation (19) since the first estimate is 
closer to the exact solution. The process is, in fact, similar to the 
iterative use of the Milne -Eddington approximation ([2], p. 87) to 
achieve a solution of the Milne equation. 

Continued iteration "beyond this stage is more difficult. A simpler 
process is available, however, in which integration of products of the 
transcendental functions is avoided while the accuracy is increased. 

To this end, in view of the known antisymmetry about g = 1^/2, we add 
a cubic term to the linear expression for cp(g). Thus, we put 

<p(6) = A (i - t0 + 7 (l - 70 (46) 


where the slope A at g = g^/2 is given by (45a). A first approxima- 
tion is found by inserting the linear form (45b) in the integral equa- 
tion (26) to find 



= A 0 ~y)+^ (1- Ag L )[U£(g L -g) -E 2 (g)] + | [E s (g) - E 3 (g L - g) ] 

(47) 

A value for 7 is determined by equating expressions (^6) and (V7)j 
evaluated at | = 0; 

.7 ^ = - Q0) 9^(0) + | Ag L 
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Now pat 


<P^(fi) = A (t - —) + (i - PT 


in the integral equation (26) and find 


cp( £ )(|) = i [E 2 (| L - - E 2 (|)3 + | [ 9 (2) (l 1 )E 1 (U - 


iJHii 

(48a) 


There results from this integration 


cp( 2 )(|) = qp (l) U) + | r (l) H(|^ L ) 


(48h) 


where 


h(s,s l ) = 


5a. - -f) Eidi - iiDdi. 


3 3 

“ 2 ( 6 " t) + 4 6 “ t) + (t) CE£(i) ■ Es(| l ■ l)] 


[e 3 (i) -e 3 (i l - i)] + 3 I l [e 4 (i) - e 4 (i l - i)] 


+ 6[e 5 (i) - EsUt. - i)] 


This process can he repeated by setting up the scheme 


■ ■■ ■IHMIHII ai mm mi mi in* inn mini i i n iii i ii 


? (3) (!) = 


r (2) (i - l 4 




,(2) = 


-CfJ 


,3 r 


cp (2) (0) + | A| l 


cp (3) (l) = cp (l) (i) + | 7 (2) H(|,! l ) 


cp (n) (i) 


,(n-i) _ 






-(^) cp (n ' l} (0) + | A| l 


> ( 49 ) 




and increasing n -until successive results agree to the desired number 
of figures. 

When a solution for cp(g) has been obtained having the required 
accuracy, the quantities q/ (q£ - q^), K 0 , K x can be calculated accord- 
ing to their definitions in equations ( 27 ) and (32). For this purpose, 
it was thought to be sufficiently accurate to use the cubic expression 


9 (n) (!) =*(l - + 7 (n - l) (i 

and evaluate the integrals analytically, rather than to perform a numeri- 
cal integration using Cf/ n ^(f-). The results are, using superscripts to 
indicate correspondence with the iterates shown in equations (49 ) ' 
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= A | 6 l [Es( 0) + Ea(6 L )] + [E 4 (0) - E 4 (| L )jj 
+ A (n_l) {- (^) 3 [E 3 (0) + E 3 (| l )] + 3 (^) 2 [E 4 (0) - E 4 (| l )] 



[E 5 (0) + Es( 6 l )] + 6 [E s ( 0) - Ee(g 



( 50 a) 


K! (ri) U L ) = A E 3 (^) - [e 4 (0) - E 4 } 

+ 6 (j 1 ) - 6 Eq( 0) - E 6 | (50b) 

Once Kq or K x is known, the flux follows from (32a) or (32h). 
Finally, numerical values for the quantity (q£ - q£) can he found hy 
using the expression (50a) for Kq for the calculation of A accord- 
ing to ( 29 ) and then using equation ( 31 a). 

IRESENTATION AND DISCUSSION OF RESULTS 

In this section the analytic solutions previously derived will be 
shown graphically and some formulas of special interest will be given. 
Attention is directed first to flux since reasonable accuracy is not 
difficult to achieve. 
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The ratio q/(q£ - Q^) &oes not distinguish betveen black -vail 
emission and the more general case. Hence, the two formulas, equa- 
tions (40a) and (45d), can be compared directly vith the numerical 
solution for black vails in [4] and [8] and also vith a recent exact 
calculation of Hottel [ 13 ]. Figure 2 shovs our predictions. The solid 
curve gives results from the iterative use of the cubic approximation 
to cp(g) as outlined folloving equation (^6). Through the range 
0 < < 2 there is excellent agreement vith published results and the 

solid curve appears to be a valid criterion to judge the merit of our 
other approximations. Unfortunately, discernible differences from the 
references can be detected in the range 3 < ^ but comparable dif- 

ferences also exist betveen the references themselves. Since there is 
complete agreement at g^ = 2 and little doubt as to the approximation 
at |^ = 10 , the slight deviations in the intermediate values are 
possibly attributable to the draftmanship of the published curves. 

The practical consequences are slight, hovever, and the present results 
appear valid to about tvo digit accuracy. The short -dash curve is 
equation (40a) (for n = 3/2) and is the interpolation formula proposed 
by Probstein [6]. Equation (45d) supplies another interpolation formula 
vhich affords a better fit on the lover end of the g-^ scale but at a 
sacrifice of algebraic simplicity. In the range 2 < g^ < 6 the inter- 
polation formulas yield about the same accuracy but differ in the sign 
of the error. For g-^ » 1 all results coalesce. 
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The next objective is to relate the magnitude of q£ - to the 
wall boundary conditions for the general case of nonblack walls. Equa- 
tions (40b) and (45e) supply this information for the two linear forms 
of the emission function. Considerable simplification is achieved, 
however, if the walls are opaque or if they are alike in material prop- 
erties. First, we shall restrict attention to opaque walls where 
1 - r Q = e Q , 1 " r L = e L* From equations (29), (31a) , and (32) one gets 


% ~ _ 1 

oT l 4 - 0T O 4 1 + [-q/(q£ - qj)][(l/e 0 ) + (l/e L ) - 2] 


(51) 


Substitution into the denominator of the right member, from equa- 
tions (40a) or (J+5&), provides the desired expression. For example, 
equation (40a) with n = 3/2 yields 

q L ~ = 1 + 3i L A __ ^ 

oT l 4 - oT 0 4 3I L A + [(1/eo) + (l/e L ) - 1] 

For large and small ^ equation (52) reduces to more familiar 
relations . 

Figures 3 show, for opaque walls, the function q/(oTj j 4 - crT 0 4 ) 
as given by the different approximations. The solid curves are pre- 
sumably exact to the accuracy with which they can be read. The short - 
dash curves axe, in effect, generalizations of Probstein* s interpolation 
formula and, from equations (ijOa) and (52), are given by the relation 


-I = 1 

oT L 4 " 3! l A + [(l/e 0 ) + (l/e L ) - 1] 


(53) 
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The remaining curves (long dash) result from use of expression (45d) 
in equation ( 51 )* The dashed curves are not shown in all cases; dif- 


ferences from the solid curve can he estimated from the results shown, 
and the fact that the three approximations are indistinguishable for 
the lowest curve on each graph. 

The relative merits of the simplest approximation, that is, the one 
associated with the exponential kernel, appear yiow to have been fully 
established insofar as predictions of flux are concerned. We therefore 
write, finally, the most general formula based upon this approach: 


6 o oT o 


4 


e L^L _ 

L(1 - r Ti ) (l - r 0 )J 


= 1 / 


L(1 - r L ) (1 - r Q ) 


- 1 


(54) 


Equation (54) embraces and extends the previous interpolation formulas. 

The actual distribution of emission, or of cp(|), provides the most 
exacting demands upon the approximations under lying the predictions. In 
figure 4 the curves were calculated by two methods. The solid lines 
again correspond to repeated iteration of the cubic expression; the 
other curves are given by equation (47)* the first iteration of the 
linear relation. No attempt was made to include the linear distribu- 
tions of cp(£) predicted by the simpler methods. To the scale of these 
figures both linear distributions conform to the straight lines tangent 
to cp(g) at | The standard of comparison comes from the numer- 

ical results plotted in [4] and [8]. It does not seem possible, because 
of inconsistencies in the published figures, to draw unequivocal con- 
clusions about the accuracy. Qualitatively, the results may certainly 
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be termed satisfactory tut small local deviations exist in each. case. 
These discrepancies are at most of the order of 1 or 2 percent. 
Diminishing preciseness inevitably occurs in the vicinity of the vails. 

In figure 5 the dependence of the functions Kq, Ki, and A on 
is shovn according to the test approximations of this paper. These 
values are of use in further specific calculations of cp(g) and in the 
determination of the level of the emission function. 

In the initial part of this paper the analysis is cast in a general 
form and the problem of radiative transport betveen vails vith rather 
arbitrary characteristics is related to the calculations involving black 
vails. In the latter part of the paper the philosophy is more one of 
pragmatism; namely, that to the accuracy of the physical idealization 
and to the accuracy of many engineering estimates, it is preferable to 
retain some analytic control over the results. In particular, the use 
of the simplifying assumptions leading to the use of the exponential 
kernel yields simple formulas of practical interest. Analytic methods 
are shown also to produce solutions of reasonable accuracy. The mathe- 
matical nature of the problem of radiative transport is of considerable 
interest to theorists, and, in the ultimate degree of perfection, numer- 
ical calculations appear to be required. Such methods are essential in 
fixing a standard of excellence and need to be understood if restric- 
tions are to be relaxed and estimates of the effects of frequency- 
dependent absorptivity and emissivity are to be studied. On the other 
hand, a good grasp of the qualitative character of the solution is 


3 ^ 



important in the extension to less idealized cases and algebraic 
formulas may be preferable to purely numerical results in filling this 
need . 

Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field , Calif., Sept. 3, 1963 
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Figure 2.- Dimensionless flux, as calculated by three methods of 
approximation, showing dependence on optical thickness . 








Figure 4.- Approximations of universal function <p( i/ Ip) for different 

optical thicknesses. 
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